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Abstract 

Predicting behavior of large-scale biochemical metabolic networks represents one of the 
greatest challenges of bioinformatics and computational biology. Approaches, such as flux 
balance analysis (FBA), that account for the known stoichiometry of the reaction network 
while avoiding implementation of detailed reaction kinetics are perhaps the most promising 
tools for the analysis of large complex networks. As a step towards building a complete 
theory of biochemical circuit analysis, we introduce energy balance analysis (EBA), which 
compliments the FBA approach by introducing fundamental constraints based on the First and 
second laws of thermodynamics. Fluxes obtained with EBA are thermodynamically feasible 
and provide valuable insight into the activation and suppression of biochemical pathways. 


Conservation principles impose constraints on the fluxes and chemical potentials associated with 
biochemical network reactions that are analogous to KirchofF s current and voltage laws for 
electrical networks (1). Flux balance analysis (2-13) invokes mass conservation, but does not 
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consider the chemical potential of non-equilibrium steady state chemical fluxes. The 
thermodynamic analysis presented here — energy balance analysis (EBA) — provides additional 
constraints on the system that are analogous to the voltage loop law. In addition to predicting 
network fluxes that are thermodynamically feasible, EBA' facilitates a detailed analysis of 
regulation of metabolic networks that is not available from FBA alone. 

The central flux balance conservation statement is given by the equation: 

Sf = 0, (l) 

where /e 9\ n is the vector of n fluxes occurring in the network, Se is the stoichiometric 
matrix, and m is the number of reactants in the system. The matrix S stores the stoichiometric 
coefficients associated with each flux in the network. In the above formulation both internal 
fluxes and boundary fluxes, which transport material into or out of the system, are included in 5. 
Typically, a number of inequalities are introduced to constrain the boundary fluxes depending 
upon the external media (7-10, 12, 13): 

<*,</,< A- (2) 

As implemented by Palsson and colleagues (2, 7-13), biological networks are assumed to 
optimize a certain biologically meaningful objective function, which is a linear combination of 
the fluxes: 

Z = £c,/,=c7. (3) 

1=1 

For example, in the analysis of E. coli central metabolism the objective is constructed from the 
production of biosynthetic precursors required to generate biomass (8, 10). Another example 
uses maximization of ATP production to simulate mitochondrial energy metabolism (11). 
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Regardless of the application, optimization of a linear objective function (Eq. 3). together with 
the equality constraints (Eq. 1) and the inequality constraints (Eq. 2). represents a linear 
optimization problem, which can be solved with linear programming (14). 


The power of FBA is illustrated by the tour de force assembly of the complete stoichiometry of 
the known reactions in E. coli central metabolism provided by Edwards and Palsson (8, 10). 
Edwards and Palsson show that the flux balance simulation of the organism-scale metabolic 
network predicts the metabolic capabilities of E. coli (8, 10). Under various external constraints 
(e.g. aerobic vs. anaerobic growth) FBA can distinguish between genes essential and not 
essential for growth in 68 of 79 cases studied (8). This predictive capability of FBA is striking 
considering the tremendously complex problems that can be modeled using little or no free 


parameters. 


Lacking from FBA is the explicit consideration of the energy balance and thermodynamics of the 
network reactions. Since biochemical networks are comprised of multiple-species reactions, 
energy balance loop equations cannot be obtained from topological loops in the network, as is 
done in electrical circuit analysis (1). However, we will show that the energy balance equations 
are obtained from an analysis of the network stoichiometry. 


If we combine redundant fluxes and remove the columns from S that correspond to boundary 
fluxes, the remaining matrix, denoted S', represents the complete set of possible internal 
chemical transformations. Using the singular value decomposition (14), S' can be decomposed 
as: 
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S' = AAB T , 


(4) 


where A has the following form: 

'4 ••• 0 0 ••• 0 " 

a= : **. i i : . (5) 

o - 4 o ••• o 

The first m columns of A form a diagonal matrix where the diagonal entries are the singular 
values A , . The entries of the remaining columns are all equal to zero. If n c is the number of 

columns of S' and r is the rank of S', then columns r + 1 through n c of the matrix B provide the 
(n c - r )-dimensional null space of S'. We introduce the matrix K, which stores the null space 
vectors of S' : 


l l 


l.r+! 


K = 


B, 


\.n c 


Bm,r + 1 
B__ 


( 6 ) 


and we denote the i row of K as k i . Summing the n c chemical reactions, each scaled by the 


corresponding entry of k j , results in a perfectly balanced reaction equation (with same reactants 

in equal proportion on either side of the equation). An example of this analysis for a simple five- 
metabolite network with four reactants is illustrated in Fig. 1 . 


To study network thermodynamics, we consider the vector A// that lists the n c chemical potential 

differences associated with the reaction fluxes. Since each k t provides weights for exactly 
balancing the chemical reaction equations (see Fig. 1), solutions to the equation 

*A// = 0 (7) 
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balance the global free energy of the network. Eq. 7 is a statement of conservation of energy, and 
hence follows from the first law of thermodynamics (15). The second law of thermodynamics 
takes the form of an inequality constraint for each flux, 

0 . ( 8 ) 

which ensures that each reaction dissipates energy. Note that the first law imposes an equality 
while the second law imposes an inequality, as expected. Eqs. 7 and 8 represent thermodynamic 
constraints that should be considered in addition to the flux balance constraints. Eq. 8 provides a 
link that couples mass balance and energy balance, and constrains the feasible flux space more 
strictly than Eqs. 1 and 2 alone. 

We have obtained the stoichiometric matrix provided by Edwards and Palsson (10), used to 
represent the flux balances in the E. coli central metabolism. The web-posted supplementary 
information (10) provides detailed descriptions of the reaction network, which contains 953 
fluxes (735 internal; 218 boundary) acting on 536 metabolites. Using the MATLAB (The 
Mathworks Inc., Natick, MA) optimization package we reproduced the linear programming 
analysis presented in Ref. 10, and optimized the production of biomass with glucose and oxygen 
uptake constrained to be less than or equal to 10 and 15 mmol g-DW 1 h* 1 , respectively. The 
resulting flux produces a growth rate of 0.81 h' 1 on glucose minimal media. 

To compute the thermodynamic properties of this large-scale network we first combine 
redundant fluxes (e.g., the phosphofructokinase A and B reactions are combined into a single 
column of S ), resulting in n r — 617 internal reactions operating on 536 metabolites. The growth 
is the optimized under the flux balance constraints (Eqs. 1 and 2) and the constraint that the 
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energy balance equations (Eqs. 7 and 8) are satisfied. We impose the additional constraint that 
A fj. must be Finite for all nonzero fluxes. The free energy can go to zero only if the flux is zero, 
implying that a given reaction is in chemical equilibrium. This analysis predicts the same optimal 
growth rate (0.81 IT 1 on glucose minimal media) as that reported above for FBA. Yet the FBA 
prediction does not represent a unique solution to the optimization problem because 
redundancies in the metabolic network allow the optimal growth rate to occur for an infinite 
number of possible internal flux distributions. Since the EBA-constrained solution is (at least in 
this case) able to achieve the same optimal growth rate as that obtained by considering only the 
FBA constraints, the fluxes obtained by EBA represent one particular optimal solution to the 
FBA linear programming problem. However, optimal flux distributions obtained by considering 
flux balance alone are not necessarily thermodynamically feasible. The EBA constraints further 
restrict the set of feasible fluxes, and provide a more physically realistic flux distribution. 
Selected fluxes from glycolysis, TCA cycle, and respiration are tabulated in Fig. 2. Fluxes from 
the wild type on glucose media under aerobic and anaerobic conditions are labeled WT (oxygen) 
and WT (anaerobic), respectively. 

EBA further allows quantitative estimation of the chemical potentials in the system. First, we 
introduce the flux resistances, defined as r t = -A/v, / / . To avoid the unphysical situation in 

which A/v, =0 for a finite f , which is equivalent to setting the flux resistance equal to zero, we 
assume that there exists a minimum flux resistance, (which is equivalent to saying that there 
exists an upper limit to the effective reaction rate constants). Thus the second law constraint is 
modified: 

/, > 0 
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&M,Z+r min f,, f t < 0 . 


(9) 


Realistically, each flux may have a different r min . However for our purposes we find that a single 
value, = 400 kcal mol’ 2 g-DW h, produces reasonable behavior and can be used to describe 
the entire network (16). The energy balance constraint can alternatively be written in terms of the 
chemical potentials (Eq. 7) or the flux resistances: 




1 

o 

1 

n 

1 

o 

n 

1 

i 

^ . 
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= Kdiag(f)r = 0. 


( 10 ) 


With the fluxes fixed using the values obtained from the above analysis, we use quadratic 
programming to find an optimal A// that minimizes the norm of the free energy vector |A/z{ . In 
addition to the fluxes, Fig. 2 lists the predicted chemical potentials and the conductances, 
c. = r ~ l , of each reaction. The choice of = 400 kcal mol’ 2 g-DW h results in reasonable 
predictions of chemical potential. For example, A//, = -9.55 kcal mol* 1 for ATP hydrolysis. 

The reaction conductance provides a measure of the activation level of the pathway. If c ( =0 

then the associated enzyme(s) is not present or is deactivated. Increases in flux or conductance, 
at a fixed free energy, indicate that a pathway is up regulated, at either the expression level or the 
post-translational level. By changing the boundary constraints, we can study how the metabolic 
network responds to changes in substrate. In Fig. 2 we compare the predicted EBA fluxes and 
free energies for the WT cell grown under aerobic and anaerobic conditions. We identify a 
pathway as “down regulated” (colored blue) if the flux conductance decreases by a factor of 4 or 
more, and “up regulated” (colored green) when the conductance increases by more than a factor 
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of 4. Based on this analysis we identify three enzymes that require activation upon moving from 
aerobic to anaerobic conditions: fumarate reductace, pyruvate formate lysase, and pyridine 
nucleotide transhydrogenase. Other pathways show increases or decreases in flux compared to 
aerobic growth. However these differences can be accounted for by changes in the free energy 
and thus do not necessarily require regulation. 

Following the work of Edwards and Palsson (8, 10) we studied the effects of gene knockouts on 
the metabolic capabilities of E. coli. We found that zwf pgi, and end knockouts result in the 
same predicted phenotype (Fig. 2), with only two up regulations compared to WT: succinate 
dehydrogenase and pyridine nucleotide transhydrogenase. Again, a number of predicted fluxes 
that differ from the WT do not require up regulation of the associated enzymes. The situation is 
similar for pyk and pgi knockouts. The pyk knockout requires an up regulation of 
phosphoenolpyruvate carboxylase to maintain growth of almost 99% that predicted for WT. No 
other significant regulations are predicted. The pgi knockout analysis predicts one significant up 
regulation and three down regulations and a similar growth rate. Thus, we find that much of the 
capacity for metabolic control is built in to the wild type expression and activation levels. Few 
regulatory steps are required when nonessential genes are knocked out; the network is robust and 
tolerant to errors and deletions (17,18). 

However, the situation is different when essential genes are knocked. The eno and pfk genes are 
examples of genes that FBA falsely identifies as nonessential to growth on glucose minimal 
media (8). The fact that eno and pfk are essential for growth in vivo (19) may be related to the 
demands that the knockouts of these genes place on metabolic regulatory mechanisms. These 
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demands are much heavier than those imposed by nonessential knockouts. Our analysis predicts 
that maintaining growth with these knockouts requires a greater number of pathway regulations 
than for the nonessential knockouts. Of the 64 reactions listed in Fig. 2, 15 are predicted to be 
upregulated and 15 downregulated for the pfk knockout. Thus, the predicted phenotype is 
markedly different from the wild type for nearly half of the reactions associated with glycolysis, 
TCA cycle, and respiration, while the pyk and pgi knockouts differ very little from WT. 

These predictions are based on one major simplification — that the entire network can be 
characterized by a single or equivalently a maximum pathway conductance. More 

realistically, constraints could be placed on each pathway based on a priori knowledge of the 
biochemistry. One promising aspect of EBA is that it is possible to incorporate as much, or as 
little, as is known about the individual reactions. For example, if we know that the ratio 
[ATP]/[ADP] in the cell is greater than the equilibrium ratio (20), then we know that the free 
energy of ATP hydrolysis satisfies the inequality: 

&Matp->adp — ATP-tADp = — 7.3 kcal mol 1 . (1 1) 

In general, consider the reaction A— > B . If the ratio [A]/[B] is greater than 1, then the free 
energy of that reaction is constrained A//, < A p." . Alternatively, if the concentration ratio is less 
than 1, then A /u t > A p." . If the reactant concentrations are known, then the constraint becomes an 
equality: = -A: s rin(Af r(7 [A]/[B]), where k B is Boltzmann’s constant, T is the absolute 

temperature, and K tn is the equilibrium constant of the reaction. 
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Together, flux balance and energy balance provide basic laws for the analysis of biochemical 
networks. These laws, akin to the basic principles of circuit theory for electrical networks, make 
the analysis and design of large-scale biochemical systems practical. The engineering approach 
to analysis and design of such complex systems is the identification of modular components that 
are separable within the system (21). Flux balance and energy balance circuit theory provides a 
basis for understanding how these modules function and interact. 
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Figure Legends 


Figure 1: Illustration of energy balance equations for a network of Five reactions. The first step is 
determination of the stoichiometric matrix from the reactions in the network. For this example, 
the stoichiometric matrix has rank r = 3 , resulting in 2 independent null space vectors. The null 
space vectors provide mutually orthogonal solutions to S'f = 0. In addition, the null space 
vectors are the rows of the matrix K, the energy balance matrix for which K AG = 0. 

Figure 2: Regulation of reactions in glycolysis, TCA cycle, and respiration. Predicted 
biochemical fluxes (in mmol g-DW' 1 h' 1 ), chemical potentials (in kcal mol' 1 ), and conductances 
(in ??) are reported for 64 central metabolic reactions, for wild type (WT) aerobic and anaerobic 
growth, and for zwf, pgl, gnd , pyk, pgi, eno, and pfk knockouts. Predicted growth rate for each 
case is reported in units of hr *. Green and blue shading indicate up and down regulation relative 
to WT, respectively. Yellow indicates that a reversible flux has changed direction. Gray shading 
indicates that the flux cannot be identified as a site of regulation because, although the 
conductance has changed moderately (by less than a factor of 10), the magnitudes of the free 
energy and the flux have both either increased or decreased. .“N.P.” denotes enzymes that are not 
present under given conditions. “Eq.” refers to reactions that are at equilibrium, for which the 
conductance cannot be estimated. All simulations are performed for glucose minimal media. 
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Stoichiometric Matrix: 



Stoichiometric Balance Equations: 

AT, , • (rxn 1 ) + AT, 2 • (rxn 2) + AT, 3 • (rxn 3) + AT, 4 • (rxn 4) + AT, 5 • (rxn 5) 


= 2A + 2B + C + D 2A + 2B + C + D 

K 2 , • (rxn 1 ) + A r 2 2 • (rxn 2) + AT 3 3 • (rxn 3) + 4 • (rxn 4) + AT 5 5 • (rxn 5) 

= 2A + 2B + C + D^>2A + 2B + C + D 
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0 

0 
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o 

0 

0 

Ea 
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1.5699 

■ 


■ 


m 

■ 
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25 
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25 
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0 
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0 

0 
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25 
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25 

9.6646 
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25 
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-4.3601 

25 

1? 003 

-7.1612 

25 

9.6546 

-3.0418 

25 
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-4.3601 

25 

17.903 
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26 

3.9011 

-1.5605 

25 

4.0624 

-1.633 

2.5 

8.3737 

-25405 

25 

5.0349 

-23730 

25 

6 8042 

-27217 

2.5 

11.517 

-4.6069 

25 

3.1375 

-1.256 

25 

3 526 

-1.4104 

25 

58706 

-23482 

25 

27074 

-1.110 

25 

3.2782 

-1.3113 

25 

5.6465 

-2.2586 

25 

3.1263 

-1.2606 

2.5 

3.5170 

-1.4072 

2A 

5.8633 
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25 

0.8665 

-5,3201 

0.1629 

■ 


■ 


■ 
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-6.3201 

01629 







0 

-0.0635 

0 
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0 

0 
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0 
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0 

0 
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0 
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0 

0 
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-03456 

25 

0.6294 

-0.2518 
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0.8638 
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25 
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-0.2516 
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0 

-16.532 

0 

0 
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0 
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0 

0 
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0 

0 
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0 

0 

Ea 

0 

0 
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0 

1.0607 

•0.4239 

25 

4.3005 

-1.7238 

25 

0 

0 

Ea 

4,7062 

-1.6821 

25 

7.5525 

-3.021 

2.5 

-02601 

3.2249 


4.7052 

•1.6821 

2.5 

7.5525 

-3.021 

25 

-02601 

02249 
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0 

0 

25 

0 

0 
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0 

0 

0 

Ea 

■i 


■ 

22606 
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0 
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0 

-6.9357 

0 
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-3.3935 

0 

0 
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0 

0 

-6.3350 

0 

0 
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0 

0 
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0 

0 

•6 3350 

0 
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29.751 -11.9 


05456 

-25.060 0.0218 

-00999 
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u 

Z W PW l Wj : 

0383 

-16.313 

0.0209 


0.2572 

-0.1029 2.5 
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-0.0362 

25 
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02534 

-0.1014 

25 
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25 

0 
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0 

0 

Ea 

o 

o 

m 

•P 
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-0.7564 

25 

7.281 

-2.9124 

25 

0 

-14.325 0 


0 

-6 5628 

0 

0 

•6.1576 

0 

-44.049 

28.649 1.5375 

-1.1514 

04605 

25 

-37.74 36.660 1.0202 

■44.509 

16.679 

23629 

-45.35 

16.14 

25 

7.6 

-0.5407 0.7956 

7 8 

-11.111 

0664 

7 6 -12.223 06216 

7.6 

-9.2262 

1.2207 

7.6 

•6.0467 

1.2560 

0.31 

-0.32 2.5 

0.285 

-0.11 

2.5 

0.734 -0.31 2.5 

0.798 

-0.32 

2JS 

0.798 

-0.32 

2.5 


0 

-72.159 

0 

29.819 

•11.928 

25 

0 

-30.307 

0 

0 

-40.382 

0 

30 

-12 

2.5 

0 

-20.604 

0 


29 632 -11.933 
0 -23.704 


29.995 -11.996 


2.5 01641 -0.0737 2.5 

25 0 0 Eq. 

0 09271 -25.913 00359 


2.5 01665 -0.0666 25 

Eq. 22.204 -6.6619 2.5 


7.6 -17.209 04416 


>.264 16.109 25 

7 6 -6 4622 1.1761 
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I Knock out 


15 




